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The detection of quasi-periodic sources of gravitational waves requires the accumulation of signal- 
to-noise over long observation times. This represents the most difficult data analysis problem facing 
experimenters with detectors like those at LIGO. If not removed, Earth-motion induced Doppler 
modulations and intrinsic variations of the gravitational-wave frequency make the signals impossible 
to detect. These effects can be corrected (removed) using a parameterized model for the frequency 

QO ■ evolution. In a previous paper, we introduced such a model and computed the number of independent 

0^ ' parameter space points for which corrections must be applied to the data stream in a coherent search. 

Q\ ■ Since this number increases with the observation time, the sensitivity of a search for continuous 

gravitational-wave signals is computationally bound when data analysis proceeds at a similar rate 
' to data acquisition. In this paper, we extend the formalism developed by Brady et al. [Phys. 

Rev. D 57, 2101 (1998)], and we compute the number of independent corrections N P (AT, N) 

^"*^ ' required for incoherent search strategies. These strategies rely on the method of stacked power 

spectra — a demodulated time series is divided into N segments of length AT, each segment is 
Fourier transformed, a power spectrum is computed, and the N spectra are summed up. This 
I ' method is incoherent; phase information is lost from segment to segment. Nevertheless, power from 

a signal with fixed frequency (in the corrected time series) is accumulated in a single frequency 
■ bin, and amplitude signal-to-noise accumulates as ~ N 1 ^ 4 (assuming the segment length AT is held 

t—{ ' fixed). We estimate that the sensitivity of an all-sky search that uses incoherent stacks is a factor of 

2-4 better than achieved using coherent Fourier transforms; incoherent methods are computationally 

O^l ' efficient at exploring large parameter spaces. We also consider a two-stage hierarchical search in 

which candidate events from a search using short data segments are followed up in a search using 
longer data segments. This hierarchical strategy yields another 20-60% improvement in sensitivity 
in all-sky (or directed) searches for old (> lOOOyr) slow (< 200Hz) pulsars, and for young (> 40yr) 
fast (< 1000Hz) pulsars. Assuming enhanced LIGO detectors (LIGO-II) and 10 12 flops of effective 
computing power, we also examine the sensitivity to sources in three specialized classes. A limited 
>' j ' area search for pulsars in the Galactic core would detect objects with gravitational ellipticities of 

(3JT), e>5x 10 -6 at 200 Hz; such limits provide information about the strength of the crust in neutron 

stars. Gravitational waves emitted by the unstable r-modes of newborn neutron stars would be 

. ^ | detected out to distances of ~ 8 Mpc, if the r-modes saturate at a dimensionless amplitude of order 

unity and an optical supernova provides the position of the source on the sky. In searches targeting 
!_l ' low-mass x-ray binary systems (in which accretion-driven spin up is balanced by gravitational-wave 

spin down), it is important to use information from electromagnetic observations to determine 
the orbital parameters as accurately as possible. An estimate of the difficulty of these searches 
suggests that objects with x-ray fluxes which exceed 2 x 10 _8 erg cm~ 2 s _1 would be detected using 
the enhanced interferometers in their broadband configuration. This puts Sco X-l on the verge of 
detectability in a broad-band search; in this case amplitude signal-to-noise might be increased by 
~ 5-10 by operating the interferometer in a signal-recycled, narrow-band configuration. Further 
work is needed to determine the optimal search strategy when we have limited information about 
the frequency evolution of a source in a targetted search. 

PACS numbers: 03.70.+k, 98.80.Cq 



I. INTRODUCTION 

The detection of periodic sources of gravitational waves 
using the LIGO, or similar, gravitational-wave detectors, 
is seemingly the most straightforward data analysis prob- 
lem facing gravitational wave astronomers. It is also the 
most computationally intensive; extremely long observa- 



tion times will be required to have any chance of detecting 
these signals. Searches for periodic (or quasi-periodic) 
sources will be limited primarily by the computational 
resources available for data analysis, rather than the du- 
ration of the signals or the lifetime of the instrument. 
For this reason, it is of paramount importance to explore 
different search strategies and to determine the optimal 
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approach before the detectors go on line at the end of 
the century. In a previous paper JjJ , hereafter referred to 
as Paper I, we presented a detailed discussion of issues 
which arise when one searches for these sources in the de- 
tector output. Using a parameterized model for the ex- 
pected gravitational wave signal, we presented a method 
to determine the number of independent parameter val- 
ues which must be sampled in a search using coherent 
Fourier transforms (which accumulate the signal to noise 
in an optimal fashion) . The results were presented in the 
context of single-sky-position directed searches, and all- 
sky searches, although the method outlined in Paper I 
is applicable to any search over a specified region of pa- 
rameter space. Livas 0, Jones (3) and Niebauer et al. ||] 
have implemented variants of the coherent search tech- 
nique without the benefit of the optimization advocated 
in Paper I. 

In this paper, we discuss alternative search algo- 
rithms which can better detect quasi-periodic gravita- 
tional waves using broadband detectors. These algo- 
rithms achieve better sensitivities than a coherent search 
with equivalent available computational resources. This 
improvement is accomplished by combining coherent 
Fourier transforms with incoherent addition of power 
spectra, and by using hierarchical searches which follow 
up the candidate detections from a first pass search. 

The most likely sources of quasi-periodic gravitational 
waves in the frequency bands of terrestrial interferometric 
detectors are rapidly rotating neutron stars. We use these 
objects as guides when choosing the scope of the exam- 
ple searches considered below. Nevertheless, the search 
algorithms are sufficient to detect all sources of continu- 
ous gravitational wave signals provided the frequency is 
slowly changing. 

A rotating neutron star will radiate gravitational waves 
if its mass distribution (or mass-current distribution) is 
not symmetric about its rotation axis. Several mech- 
anisms which may produce non-axisymmetric deforma- 
tions of a neutron star, and hence lead to gravita- 
tional wave generation, have been discussed in the lit- 
erature |p|-|lO| ■ A neutron star with non-zero quadrupole 
moment which rotates about a principle axis produces 
gravitational waves at a frequency equal to twice its ro- 
tation frequency. Equally strong gravitational waves can 
be emitted at other frequencies when the rotation axis 
is not aligned with a principal axis of the source |flp|- 
If the star precesses, the gravitational waves will be pro- 
duced at three frequencies: the rotation frequency, and 
the rotation frequency plus and minus the precession fre- 
quency (§. 

For concreteness, we consider a model gravitational- 
wave signal with one spectral component. This is not a 
limitation of our analysis since the search strategy pre- 
sented below is inherently broadband; it can be used to 
detect sources which emit gravitational waves at any fre- 
quency in the detector pass-band. Additional knowledge 
of the spectral characteristics of a signal might allow us to 
improve our sensitivity in the case when multiple spectral 



components have similar signal-to-noise ratio. In such a 
circumstance, a modified search algorithm would sum the 
power from each frequency at which radiation would be 
expected. In a background of Gaussian noise, the sensi- 
tivity would improve as (number of spectral lines) 1 / 4 for 
only a moderate increase in computational cost. 

Finally, we mention several other works which con- 
sider searching for quasi-periodic signals in the output 
of gravitational wave detectors. Data from the resonant 
bar detectors around the world has been used in searches 
for periodic sources. New et al. |l2j have discussed is- 
sues in searching for gravitational waves from millisecond 
pulsars. Krolak [|l3| and Jaranowski et al. |[4|,[l5| have 
considered using matched filtering to extract information 
about the continuous wave sources from the data stream. 
Finally, work is ongoing in Potsdam to investigate line- 
tracking algorithms based on the Hough transform ffg) ; 
this technique looks quite promising, although we must 
await results on the computational cost and statistical 
behavior before we can make a detailed comparison to 
the techniques described in this paper. 



A. Gravitational waveform 

The long observation times required to detect contin- 
uous sources of gravitational waves make it necessary to 
account for changes in the wave frequency; the physi- 
cal processes responsible for these changes, and the as- 
sociated time scales were discussed in Paper I. In addi- 
tion, the detector moves with respect to the solar system 
barycenter (which we take to be approximately an iner- 
tial frame) , introducing Doppler modulations of the grav- 
itational wave frequency. To account for these two effects, 
we introduce a parameterized model for the gravitational 
wave frequency f(t; A) and phase 4>(t; A) = 2n Jf(t; A) dt 
measured at the detector: 



f(t;X) = f (l + --n) [l+J2f> 
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Here fo is the initial, intrinsic gravitational-wave fre- 
quency, x(t) is the detector position, v(t) is the detector 
velocity, n is a unit vector in the direction of the source, 
and fk are arbitrary coefficients which we call spindown 
parameters. (We refer the reader to Paper I for a detailed 
discussion of this model and its physical origin.) The vec- 
tor A denotes the search parameters — the parameters of 
the frequency model that are (generally) unknown in ad- 
vance. In the most general case that we consider below, 
the search parameters include frequency fo, the polar 
angles (9, if) used to specify n, and the spindown param- 
eters f k : 

A = (A , A 1 , A 2 , A 3 , A 4 , . . .) = (f ,e, <p, fi,h, ■■■)■ (1.3) 
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We note that the parameter A = fo defines an overall 
frequency scale, whereas the remaining parameters define 
the shape of the phase evolution. It is convenient to 
introduce the projected vector A = (A 1 , A 2 , A 3 , A 4 , . . .) of 
shape parameters alone. 

The strain measured at the interferometer is a linear 
combination of the + and x polarizations of the gravita- 
tional waves, and can be written as 



Recent work 20-23 has suggested that new-born 



h(t;X) = Re[^c- i ^ t;A )+*]] 



The time-dependent amplitude A and phase 'J depend 
on the detector response functions and the orientation of 
the source; they vary gradually over the course of a day 
(see references [[7p^|). In what follows, we treat A and 
^ as constants. Our analysis may be generalized to in- 
clude the additional phase modulation; however, this ef- 
fectively increases the dimension of the parameter space 
by one and the number of points that must be sampled 
by ~ 4, which translates into a reduction in relative sen- 
sitivity of ~ 6%. 



B. Parameter ranges 

The computational difficulty of a search for quasi- 
periodic signals depends on the range of parameter values 
that are considered in the search. The intrinsic grav- 
itational wave frequency /o ranges from (near) zero to 
some cutoff frequency /max- If gravitational waves are 
emitted at twice the rotation frequency, theoretical esti- 
mates [ [l7|]lgl l suggest that 



/max ^ 1.2 kHz to 4 kHz 



(1.5) 



depending on the equation of state adopted in the 
neutron star model. Observational evidence — the co- 
incidence of the periods of PSR 1937+21 and PSR 
1957+20 — favors the lower bound on gravitational wave 
frequency / max — 1.2 kHz [p^|. The spindown param- 
eters fj are allowed to take any value in the range 
\fj\ < (l/ T min) J where r min ~ /// is the characteristic 
time scale over which the frequency might be expected 
to change by a factor of order unity. Observations of ra- 
dio pulsars provide rough guidance about the time scales 
Tmin- In Paper I we considered two fiducial classes of 
sources which we denoted: (i) Young, fast pulsars, with 
/max = 1000 Hz and T m i n = 40 yr, and (ii) old, slow 
pulsars, with / max = 200 Hz and T nl - m = 1000 yr. To 
facilitate direct comparison with the achievable sensitiv- 
ities quoted in Paper I, we again use these two classes to 
present our results. 

The two extremes of sky area to be searched are: (i) 
zero steradians for a directed search in which we know 
the source location in advance; e.g. a supernova remnant, 
and (ii) 4-7r steradians for an all-sky search. We consider 
both of these cases, as well as the intermediate case of a 
0.004 steradian search about the galactic center. 



rapidly spinning neutron stars may evolve on a time 
scale of months rather than decades, radiating away most 
of their angular momentum in the form of gravitational 
waves within a year. These sources may be loud enough 
to be detected in other galaxies, in which case optical 
detection of a supernova can serve as a trigger for a tar- 
geted search. Therefore we consider the case of a directed 
search for sources with frequencies of / ma x=200 Hz and 



(1.4) evolution time scales of r m i n =l yr. 



A final class of sources that we consider are accret- 
ing neutron stars in binary systems. Several such binary 
systems have been identified via x-ray observations; the 
rotation frequencies of the accreting neutron stars are in- 
ferred to be ~ 250-350 Hz (/ max =700 Hz). Bildsten || 
has argued that these accreting objects in low mass x-ray 
binaries (LMXB's) may emit detectable amounts of grav- 
itational radiation. Since the positions of these sources 
are well localized on the sky by their x-ray emissions, 
the earth-motion induced Doppler modulations of the 
gravitational waves can be precisely determined. The 
difficulty with these sources is the unknown, or poorly- 
known, orbits of the neutron stars about their stellar 
companions, and the stochastic accretion- induced vari- 
ations in their spin. We have estimated the size o f these 
effects, and outlined a search algorithm in Sec. VHC. 
These issues deserve further study in an effort to improve 
the search strategy. 



C. Search Technique 

In searches for continuous gravitational waves, our sen- 
sitivity will be limited by the computational resources 
available, rather than the duration of the signal or the 
total amount of data. Therefore the computational ef- 
ficiency of a search technique is extremely important. 
For example, matched filtering (convolution of noise- 
whitened detector output with a noise-whitened tem- 
plate) may detect a signal with the greatest signal-to- 
noisc ratio for any given stretch of data; however, it be- 
comes computationally prohibitive to search over large 
parameter spaces with long data stretches. A sub- 
optimal, but more efficient, algorithm might achieve the 
best overall sensitivity for a fixed amount of computa- 
tional resources. We present two possible search strate- 
gies to accumulate signal to noise from the data stream. 

Central to both of these methods is the technique we 
adopt to demodulate the signal. We can remove the ef- 
fects of Doppler and spindown modulations by defining 
a canonical time coordinate 

fc+i 
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with respect to which the signal, defined in Eq. (|1 
perfectly sinusoidal: 
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lecord / and X 
for any bins above • 
threshold p 



Acquire data series 
of length T 



Divide data series 
into A? stacks 



Choose large 
patch X 



Correct phase in each 
stack by resampling 
according to X 



FFT each stack 



Choose small patch X' 
inside large patch 



T 



Correct residual frequency 
drift by sliding spectra 
according to X' 




FIG. 1. A flowchart representation of the stacked slide 
algorithm to search for sources of continuous gravitational 
waves. Notice that the computational cost of sampling the 
fine grid is reduced by sliding the power spectra, rather than 
re-computing an FFT for each point on the fine grid. 



h(t b [t;\]) = Ae~ 27riht »^- 
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(Remember, we treat A and \& as constant in time.) The 
introduction of the new time coordinate can be achieved 
by re-sampling the data stream at equal intervals in 
tf,. The power spectrum, computed from the Fourier 
transform of the re-sampled data, will consist of a sin- 
gle monochromatic spike, whose amplitude (relative to 
broadband noise) increases in proportion to the length of 
the data stretch. In practice the data will be sampled in 



the detector frame, so that a sample may not occur at the 
desired value of tj,. Consequently, we advocate the use 
of nearest- neighbor (stroboscopic) resampling p^j . This 
method will not substantially reduce the signal to noise 
in a search provided the detector output is sampled at a 
sufficiently high frequency. (See Appendix |b[) 

When the waveform shape parameters A are not known 
in advance, one must search over a mesh of points in 
parameter space. The result of a phase correction and 
Fourier transform will be sufficiently monochromatic only 
if the true signal parameters lie close enough to the one of 
the mesh points. In Sec. || we rigorously define what is 
meant by "close enough" , and show how to determine 
the number of points for which corrections should be 
applied. We note that the approach of resampling fol- 
lowed by a Fourier transform has the benefit that a sin- 
gle Fourier transform automatically searches over all fre- 
quencies /o, leaving only the shape parameters A to be 
searched explicitly. Other demodulation techniques, such 
as matched filtering, must apply separate corrections for 
each value of fo in addition to the A. This increases the 
computational cost dramatically. 

A signal can also be accumulated incoherently from 
successive stretches of data by adding their power spec- 
tra |2^| . However, even if each data stretch is demodu- 
lated to sufficient precision that the power from a signal is 
focused in a single Fourier frequency bin, residual errors 
in A may cause the power to be at different frequencies 
between successive spectra. A more precise knowledge of 
the phase evolution is required to correct for this drift; 
i.e. a finer mesh in parameter space. However, once a set 
of parameter corrections AA is assumed, it is relatively 
easy to correct for the frequency drift: successive power- 
spectra are shifted in frequency by a correction factor 
A/, where A/ is computed by differencing /(i;A), in 
Eq. ( |1 . ip , between the initial and corrected guesses for 
A, as a function of the start time of each data stretch. 
Once the spectra have been corrected by A/, they can 
be added together. This accumulates signal-to-noise less 
efficiently than coherent phase corrections and FFT's, 
but is computationally cheaper. 



1. Stack-slide search 

The search techniques that we consider in this paper 
are variants on the following scheme. First, the data 
stream is divided into shorter lengths, called stacks. Each 
stack is phase corrected and FFT'ed, using a mesh of 
correction points sufficient to confine a putative signal to 
~ 1 frequency bin in each stack. The individual power 
spectra are then corrected for residual frequency drift 
using a finer parameter mesh suitable to remove phase 
modulations over the entire data stretch. The corrected 
power spectra are summed, and searched for spikes which 
exceed some specified significance threshold [p6|. The 
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complete procedure is summarized in the flowchart in 
Fig.f 




List of 
"candidate sources' 
bins / and patches % 
exceed threshold ff > 





First pass: 

Search entire range 
in / and A, using 
P (1 



threshold 




Second pass: 

Search only bins in / and 
patches in X around 
candidate sources, using 



threshold « (2) 




FIG. 2. A flowchart representation of the hierarchical 
algorithm to search for sources of continuous gravitational 
waves. It should be noted that while this approach will almost 
certainly be incorporated into the eventual search algorithm 
for gravitational waves, the real benefit of such an approach 
will be to increase the confidence in a detection made using 
some other technique. 



2. Hierarchical search 



mesh points are optimally chosen between the first and 
second passes. 
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FIG. 3. Characteristic amplitudes h c [see Eq. (3.5) in Q] 
for several postulated periodic sources, compared with sensi- 
tivities h 3 / yr of the initial, enhanced and advanced detectors 
in LIGO. (h 3 / yr corresponds to the amplitude h c of the weak- 
est source detectable with 99% confidence in iyr = 10 7 s inte- 
gration time, if the frequency and phase of the signal, as mea- 
sured at the detector, are known in advance.) Long-dashed 
lines show the expected signal strength as a function of fre- 
quency for pulsars at a distance of 8.5 kpc assuming a gravita- 
tional ellipticity e = 10~ 5 of the source (see Ref. g). Upper 
limits are plotted for the Crab and Vela pulsars, assuming 
their entire measured spindown is due to gravitational wave 
emission. The characteristic amplitude of waves from r-modes 
is also shown. These signals are not precisely periodic; rather, 
they chirp downward through a frequency band of ~ 200 Hz 
in 2 x 10 r seconds. Finally, the strength of the gravitational 
waves from LMXB's, normalized to the observed x-ray flux 
from Sco X-l, is plotted under the assumption that gravita- 
tional waves are entirely responsible for their angular momen- 
tum loss. 



We also consider a two-pass hierarchical search strat- 
egy. In this case, one performs an initial search of the 
data using a low threshold which allows for many false 
alarms. This is followed by a second pass, using longer 
stretches of data, but searching the parameter space only 
in the vicinity of the candidate detections of the first pass. 
This procedure is summarized in Fig. The advantage 
of a hierarchical search are two-fold: (i) the low thresh- 
old on the first pass allows detection of low-amplitude 
signals which would otherwise be rejected, and (ii) the 
second pass can search longer data stretches on a lim- 
ited computing budget, because of the reduced param- 
eter space being searched, thus excluding false positives 
from the first pass. For given computational resources, 
this technique achieves the best sensitivity of the strate- 
gies considered here and in Paper I, if the thresholds and 



D. Results 

T he s ensitivity 8 = 1/hth of a search is defined in 
Eq. (3.4). The threshold strain amplitude h t h is defined 
such that there is a 1% a priori probability that detector 
noise alone will produce an event during the analysis, and 
therefore is the minimum characteristic strain detectable 
in the search. We compare our results for the sensitiv- 
ity to a canonical sensitivity determined by the search 
threshold h 3/yl = 4.2yJS n (f) x 10- y Hz. This threshold 
is the characteristic amplitude of the weakest source de- 
tectable with 99% confidence in a coherent search of 10 7 
seconds of data, if the frequency and phase evolution of 
the signal are known. The relative sensitivity ro i is given 
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by Srei = hs/yr/hth', a relative sensitivity re i = 0.1 for 
a search means that a signal must have a characteristic 
amplitude h c > 10 x ft. 3 / yr to be detected in that search. 
Figure || shows h 3 / yr based on noise spectral estimates 
for three detector systems in LIGO: the initial detectors 
are expected to go on-line in the year 2000, with the 
first science run from 2002-2004; the upgrade to the en- 
hanced detectors should begin in ~2004, with subsequent 
upgrades leading to, and perhaps past, the advanced de- 
tector sensitivity. The expected amplitudes h c of several 
putative sources are also shown; we use the definition 
of h c given in Eq. (50) of Ref. §7|, and Eq. (3.5) of 
Paper I. The strengths of gravitational waves from the 
Crab and Vela radio pulsars are upper limits assuming 
all the rotational energy is lost via gravitational waves. 
The estimates of waves from the r-mode instability are 
based on Owen et al. [f23| , and those from Sco X-l are 
based on the recent analysis by Bildstcn p4j. 
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FIG. 4. A contour plot of the relative sensitivity denned 
as a function of available computing power, and 
the number of stacks in the search. The plot indicates the 
sensitivity that can be achieved using a stack-slide search for 
sources with / < 1000 Hz and r > 40 yr. The darkest shading 
represents the worst sensitivity. For fixed number of stacks 
the sensitivity improves with increasing computing power as 
expected. Notice that for fixed computing resources there is 
generally a point of optimal sensitivity; indeed the number 
of stacks at this optimal operating point should be compared 
with those given in Fig. |^. It is important to notice that 
the maximum falls off very slowly as the number of stacks 
increases. 

A reasonable long-term search strategy is that data 
analysis should proceed at roughly the same rate as data 
acquisition. Given finite computational resources and a 
desired overall false alarm probability, there is an optimal 
choice for the length of a data stretch, and the number 
of stacks that one should analyze in a given search run. 
The optimal strategy is that which maximizes the final 



sensitivity of a search subject to the constraints on com- 
putational resources and time to analyze the data. We 
have plotted the relative sensitivity of a search for young, 
fast pulsars as a function of the number of stacks and the 
available computing power in Fig. ^. The optimal num- 
ber of stacks is easily read off the plot for fixed computing 
power. Note that the maximum sensitivity in this plot 
is quite flat, especially in the regime where one is most 
computationally bound. This may be extremely relevant 
when implementing these search techniques; data man- 
agement issues may impose more severe constraints on 
the size and number of stacks than computational power 
does. This remains to be explored when the data analysis 
platforms have been chosen. 



1 



w 
S- 
o 
m 

CD 

> 
o 



(a) 



m 
f- 

w 

CD 

> 

'+■ 
CD 



0.1 



0.0 




(i) All sky 40 Yrs, 1000Hz 

(ii) All sky 1000 Yrs, 200Hz 

(iii) Directed 40 Yrs, 1000Hz 

(iv) Directed 1000 Yrs, 200Hz 



10 n 10 12 10 13 10 14 10 15 10 16 
Computing power (flops) 

1 



0.1 



0.01 




(i) All sky 40 Yrs, 1000Hz 

(ii) All sky 1000 Yrs, 200Hz 

(iii) Directed 40 Yrs, 1000Hz 

(iv) Directed 1000 Yrs, 200Hz 



10 11 10 12 10 13 10 14 10 15 10 16 
(t>) Computing power (flops) 

FIG. 5. Relative amplitude sensitivities O la \ = h 3 / yr /hth 
achievable with given computational resources, for (a) one 
pass stack-slide search strategies, and (b) two-pass hier- 
archical strategies (using the stack-slide algorithm in each 
pass). The results are presented for our fiducial classes 
of sources: (i) all-sky search for young (r > 40 yr), fast 
(/ < 1000 Hz) pulsars, (ii) all-sky search for old(r > 1000 yr), 
slow (/ < 200 Hz) pulsars, (iii) directed search for young, fast 
pulsars, and, (iv) directed search for old, slow pulsars. For 
a given computational power, we have determined the opti- 
mum observation time, and number of stacks as described in 
Sees. and 0. Thus Tith is the expected sensitivity of the 
detector for an optimal stack-slide search, with 99% confi- 
dence. 

Figure ||(a) shows the optimal sensitivities which can 
be achieved, as a function of available computing power, 
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using a stack-slide search. The results are presented for 
both fiducial classes of pulsars: old (r > 1 000 yr) slow 
(/ < 200 Hz) pulsars, and young (r > 40 yr) fast (/ < 
1 000 Hz) pulsars. In each case, we have considered both 
directed and all-sky searches for the sources. The results 
should be compared with those of Paper I, in which we 
considered coherent searches without stacking: the use 
of stacked searches gains a factor of ^2-4 in sensitivity. 

The use of a two-pass hierarchical search can further 
improve sensitivities by balancing the computational re- 
quirements between the two passes. Figure |f](b) shows 
the sensitivities achievable when each pass uses a stack- 
slide strategy. The sensitivities achieved exceed those of 
one-pass stack-slide searches by ~ 20-60%. 

The computational requirements for all-sky, all fre- 
quency surveys are sufficiently daun ting that we explore 
three restricted searches in Sec. |VIl[ (i) a directed search 
for a newborn neutron star in the young (< 1 year 
old) remnant of an extra-galactic supernova, (ii) an area 
search of the galactic core for pulsars with t > 100 yr 
and / < 500 Hz, and (iii) a directed search for an ac- 
creting neutron star in a binary system (such as Sco X- 
1). Figure || shows the relative sensitivities attainable 
in such searches. With computational resources capable 
of 1 Tflops, we expect to see galactic core pulsars with 
enhanced LIGO if they have non-axisymmetric strains of 
e > 5 x 10~ 6 at frequencies of ~ 200 Hz. Estimates of 
the characteristic strain of gravitational waves from an 
active r-mode instability in a newborn neutron star sug- 
gest that these sources will be detectable by the enhanced 
interferometers in LIGO out to distances ~ 8 Mpc; the 
rate of supernovae is ~ 0.6 per year within this distance. 
Finally, gravitational waves from accreting neutron stars 
in low-mass x-ray binary systems (LMXBs) may be de- 
tectable by enhanced interferometers in LIGO if we can 
obtain sufficient information about the binary orbit from 
electromagnetic observations. Sco X-l is on the margins 
of detectability using the enhanced LIGO interferometers 
operating in broadband configuration. We estimate that 
the amplitude signal-to-noise from these sources could be 
improved by a factor of ~ 5-10 by operating the interfer- 
ometer in a signal-recycled, narrow-band configuration. 



E. Organization of the paper 

In Sec. [n]we extend the metric formalism that was de- 
veloped in Paper I to determine the number of parameter 
space points that must be sampled in a search that ac- 
cumulates signal to noise by summing up power spectra. 
This method can then be used to compute the number of 
correction points needed in a stack-slide search. Approx- 
imate formulae, useful for estimating the computational 
cost of a search, are presented for the number of correc- 
tions needed in an all-sky search, and also in directed 
searches of a single sky position. 



Then we present the computational cost estimates, and 
determine the optimal parameters for single-pass, stack- 
slide searches in Sec. |[y] 
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FIG. 6. Panel (a) represents the relative amplitude sensi- 
tivities O rc i = /i3/ yr /ftth achievable with given computational 
resources, in three specialized searches: (i) A search for a 
new-born neutron star (whose direction is determined by ob- 
serving an optical supernova) that is spinning down by grav- 
itational wave emission via an active r-mode instability. We 
took r > 1 yr and / < 200 Hz. (ii) A search for pulsars in a re- 
gion extending 0.004 steradians about the galactic core, with 
t > 100 yr and / < 500 Hz. (iii) A source in a binary orbit, 
e.g. Sco X-l. We assume the orbit is characterized by two or- 
thogonal velocity parameters, known to within a total error of 
17(km/s) 2 ; we further assume that the frequency / < 500 Hz 
experiences a random walk typical of Eddington-rate accre- 
tion. For each of these sources, panel (b) shows h t h for initial 
(upper lines), enhanced (middle lines), and advanced (lower 
lines) interferometers in LIGO, assuming 1 Tflops of comput- 
ing power. Thus this is the characteristic amplitude of the 
weakest source that can be detected with 99% confidence us- 
ing a two pass-hierarchical search strategy. 

Section |v] presents a general discussion of hierarchi- 
cal searches for periodic sources using a single interfer- 
ometer. Schutz [ 2Sj has previously considered hierarchi- 
cal searches demonstrating their potential in searches for 
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periodic sources. The relationship between the thresh- 
old in the second stage of the search, and the thresh- 
old required in the hrst stage is discussed in detail. We 
also present the computational cost of each stage of the 
search. 



VI 



to determine 



These results are used in Sec. 
the optimal search parameters in hierarchical searches 
for our hducial classes of sources. 

Fi nally , we discuss three specialized searches in 
Sec. VII. We present a preliminary investigation of is- 
sues which arise when the gravitational wave source is 
in a binary system (e.g. an LMXB). We discuss a pa- 
rameterized model of the binary orbit, and estimate the 
number of parameter space points which must be sam- 
pled in a search for the gravitational waves from one of 
the objects in the binary. In the case where the emitter 
is accreting material from its companion, we also allow 
for stochastic changes in frequency due to fluctuations in 
the accretion rate. 

Detailed formulae for the number of points in param- 
eter space when dealing with a stacked search are pre- 
sented in Appendix |X| In Appendix |b] we discuss the 
loss in signal to noise that can occur when using nearest 
neighbor resampling to apply corrections to the detector 
output. If the data is sampled at 16 384 Hz, we demon- 
strate that this method will lose less than 1% of ampli- 
tude signal-to-noise for a signal with gravitational wave 
frequency < 1000 Hz, 



II. MISMATCH 

In a detection strategy that searches over a discrete 
mesh of points in parameter space, the search parameters 
and signal parameters will never be precisely matched. 
This mismatch will reduce the signal to noise since the 
signal will not be precisely monochromatic. It is desir- 
able to quantify this loss, and to choose the grid spacing 
so that the loss is within acceptable limits. This can be 
achieved by defining a distance measure on the parame- 
ter space based on the fractional losses in detected signal 
power due to parameter mismatch. In Paper I we derived 
such a measure in the case where the search was per- 
formed using coherent Fourier transforms; this method 
was modeled after Owen's computation of a metric on 
the parameter space of coalescing binary waveforms [ p9| . 
In this paper we extend this approach to the case of inco- 
herent searches, in which several power spectra are added 
incoherently, or stacked, and then searched for spike s. 

Let h(t; A) be a hypothetical signal given by Eq. ( |l.4| ) 
with true signal parameters A = (fo,X). If the data 
containing this signal are corrected for some nearby set 
of shape parameters A + AA, the signal will take the form 

h b (t; A, AA) = _4 e - i {2T/ot i ,[t;A]+0[t;A]-0[t;(/ o ,A+AA)]} ^ 

(2.1) 

where the subscript b is used to indicate the corrected 
waveform. In a stacked search, the data are divided into 



N segments of equal length AT, each of these segments 
is Fourier transformed, and then a total power spectrum 
Ph(f; A, AA) is computed according to the formula 



N 

E 

k=l 



ff(/;A,AA) = 2V|M/;A,AA)| 



(2.2) 



The Fourier transform of each individual segment is de- 
fined to be 



h k (f;X,AX) 



A 



rkAT 



giA0[t;A,AA] dh (2 3) 



VAT J(fc-l)AT 

where A(f> is given by 

A#; A, AA] = 2tt(/ - f )t b + <f>[t; (/„, A + AA)] - #; A] . 

(2.4) 

Here AA = (/ — /o,AA) denotes the error in match- 
ing the modulation shape parameters and the error in 
sampling the resulting power spectrum at the wrong fre- 
quency. Both of these errors lead to a reduction in the 
detected power relative to the optimum case where the 
carrier frequency and the phase modulation are precisely 
matched. 

The mismatch m(A, AA), which is the fractional re- 
duction in power due to imperfect phase correction and 
sampling at the wrong Fourier carrier frequency, is de- 
fined to be 



m(A, AA) = 1 - 



g(/;A,AA) 
H(f o ;X,0) 



(2.5) 



Remember A = (A , A) = (/, Aj, A2 . . .). Substituting the 
expressions for H from Eq. (2J2) into Eq. (|2.5| ), we find 



1 N 

(A,AA) = 1-— £M/;A,AA)| 



(2.6) 



fc=i 



It is easily shown that m(A, AA) has a local minimum of 
zero when A A = 0. We therefore expand the mismatch 
in powers of AA to find 



m(A, AA) = 5a/3(A)AA a AA' 3 + C(AA 3 ) , (2.7) 

a,f3 



where (a, f3) are summed over 0, 1, . . . , j, . . .. The quan- 
tity g a (3 is a local distance metric on the parameter space. 
This metric is explicitly given by 



g a p(X) = \d^x"d AX (im(\ AA)|aa=o , 



(2.8) 



where <9aa° denotes a partial derivative with respect to 
AA Q . It is convenient to express g a f3 as a sum of metrics 
computed for the individual stacks, that is 



/c=i 



(2.9) 
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where the individual stack metrics g^l (A) are explicitly 
given by 

(2.10) 



The phase error A0 is given in Eq. (2.4), and we use the 
notation 



Ik AT 



kAT 



(. . .)dti 



(fc-l)AT 



(2.11) 



AA=0 



In a search, we will look for spikes in the power spec- 
trum H(f; A, AA) computed from the detector output, 
that is, we will look for local maxima in the frequency 
parameter /. Therefore the relevant measure of distance 
in the space of shape parameters A is the fractional loss in 
power due to mismatched parameters AA, but after max- 
imizing over frequency. We therefore define the projected 
mismatch fi(X, AA) to be 

/i(A, AA) = minm(A, A A) = ^ 7 y(A)AA i AA J , (2.12) 
where 



7ij 



9ij 



9io9jQ 
5oo 



(2.13) 



A°=/ m 



is the mismatch metric projected onto the subspace of 
shape parameters, and / max is the maximum frequency 
that we include in the search. The meaning of the min- 
imization min f is clear from the definition of the mis- 
match in Eq. (J2.6|). 

The distance function, and in particular the metric in 
Eq. (2.13), can be used to determine the number of dis- 



crete mesh points that must be sampled in a search. Let 
V be the space of all parameter values A to be searched 
over, and define the maximal mismatch /i max to be the 
largest fractional loss of power that we are willing to tol- 
erate from a putative source with parameters in V . For 
the model waveform in Eqs. (1.1), (1-2), and ( |1.4| ) this 
parameter space is coordinatized by A = (6, (f), fx, f%, . . .) 
where 9,<fi denote location of the source on the sky, and 
fj are related to the time derivative of the intrinsic fre- 
quency of the source. Each correction point of the mesh 
is considered to be at the center of a cube with side 
2\/ Umax/n, where n is the dimension of V; this insures 
that all points in V are within a proper distance /z max of a 
discrete mesh point as measured with the metric 7^ . The 
number of patches required to fill the parameter space is 



N P (AT AT) 



U yd^~\d n \ 

(2 v /^ max /n)" 



(2.14) 



Since /t/ max is the maximum loss in detected power after 
the power spectra have been added, N P (AT, fi max , N) is 



the number of patches required to construct the fine mesh 
in the stacked search strategy described in Sec. IC. The 



coarse mesh in the stacked search strategy requires only 
that the spikes in the individual power spectra be re- 
duced by no more than /x max ; consequently, the number 
of points in such a mesh is simply N p (AT, fj, max , 1). 

We note that the average expected power loss for a 
source randomly placed within such cubical patches is 
(n) = ^m ax /3. (In Paper I we quoted an average that was 
computed for ellipsoidal patches; this is not appropriate 
to the cubical grid which will likely be used in a real 
search.) 



A. Directed search 



In most cases, the forms of Eqs. (1.2) and ( [L6| ) are 
sufficiently complicat ed t o defy analytical solution, espe- 
cially since v in Eq. (1.2) should properly be taken from 
the true ephemeris of the Earth during the period of ob- 
servation. However, for the case of a directed search, 
that is a search in just a single sky direction, the phase 
correction is polynomial in t, and the metric can be com- 
puted analytically. To a good approximation, the met- 
ric is flat — the spacing of points in parameter space 
is independent of the value of the spindown parameters 
(/it /2> • • • , fs)- For a given number s of spindow n pa - 
rameters in a search, the right hand side of Eq. ( 2.14 ) 
can be evaluated analytically. The result is expressed as 
a product Af s G s , where 



/max (AT) 



s(s+3)/2 



(u / s wv (s+1)/2 

(Umax/*) ' mm 



(2.15) 



depends on the maximum frequency / max (in Hz), the 
length of each stack AT (in seconds), the maximal mis- 
match /i max , and the minimum spindown age r m i n (in 
seconds) to be considered in the search. The dependence 
on the number of stacks N is contained in G S (N) which 
are given by: 



G (N) = 1 , 
Gi(JV) « 0.524iV , 
G 2 (N) w 0.07087V 3 , 
G 3 (N) « 0.00243JV 6 



(2.16) 
(2.17) 
(2.18) 
(2.19) 



when N 3> 4. The detailed expressions for G S (N) are 
presented in Appendix [a|. For up to 3 spindown terms 
in the search, the number of patches is then 



N P (AT N) 



max W S G S (N)} . (2.20) 
se{o, 1,2,3} 



The maximization accounts for the situation where in- 
creasing s, the dimension of the parameter space, de- 
creases the value of A/" S G S because the parameter space 
extends less than one patch width in the new spindown 
coordinate f s ; one should not search over this coordinate. 
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B. Sky search 

For signal modulations that arc more complicated than 
simple power-law frequency drift, it is impossible to com- 
pute N p analytically. In an actual search over sky po- 
sitions as well as spindown, one should properly com- 
pute the mismatch metric numerically, using the exact 
ephemeris of the Earth in computing the detector posi- 
tion. In Paper I we computed N p (Tb, ^ max , 1) numeri- 
cally, with the simplification that both the Earth's rota- 
tion and orbital motion were taken to be circular. How- 
ever, in this paper we are concerned also with the de- 
pendence of N p on the number of stacks N. This signif- 
icantly complicates the calculation of the metric and its 
determinant, and makes it necessary to adopt some ap- 
proximations in the calculation. Fortunately, the results 
of interest here are rather insensitive to errors in N p . 

In Paper I we mentioned that there are strong corre- 
lations between sky position and spindown parameters, 
thus requiring the use of the full s + 2 dimensional met- 
ric. However, these correlations are due primarily to the 
Earth's orbital motion, to which a low-order Taylor ap- 
proximation is good for times much less than a year. 
Therefore we treat the number of patches as the product 
of the number of spindown patches times the number of 
sky positions M s , computed analytically using only the 
Earth's rotational motion. We note that this approxi- 
mation is appropriate only for computing the number of 
patches; when actually demodulating the signals, the true 
orbital motion would have to be included. This approx- 
imation works well so long as the orbital residuals (the 
remaining orbital modulations after correction on this 
sky mesh) are much smaller than the spindown correc- 
tions being made at the same power in t. The residual 
orbital velocity at any power t k is roughly 



x — x (nt) k , 



(2.21) 



where is a number of order unity, M. s is the number of 
sky patches, and r = 1AU and £1 = 2n/yr are the Earth's 
orbital radius and angular velocity. When the range in 
this residual is comparable to or larger than the range in 
the corresponding spindown term fkt k , the "spindown" 
parameter space must be expanded to include the orbital 
residuals. While the range in ctk is difficult to arrive at 
analytically, we have found that assuming a maximum 
value of « 0.3 gives good agreement with the numerical 
results of Paper I (i.e. for N = 1), to within factors of 
- 2. 

One other approximation was made in computing the 
number of sky patches. We found that the measure 
ydet \^fij] for the sky position metric is almost constant 
in the azimuth <p, and has a polar angle dependence which 
is dominantly of the form sin 28. When performing the 
integral over sky positions, we approximated the measure 
by Y^deFpyiJf ~ constant x sin 29; this approximation is 
accurate to about one part in 10 4 . 



Given these approximations, the number of patches for 
a sky search is given by: 



N„ = max 

sG{0,l,2,3} 



M s jr s G s f[ ( i 



fc=0 



0.3rn k+1 T k ir 



(2.22) 



The number of sky patches M s , in the (s+2)-dimcnsional 
search, is given approximately by 



M 



2 

max 



4/i max /(s + 2) 
A = 0.014 , 

B = 0.046(AT/1 day) 2 , 
C = 0.18(AT/lday) 5 N 3 



B~ 



C~ 



-1/2 



(2.23) 

(2.24) 
(2.25) 
(2.26) 



This is a fit to the analytic result given in Appendix |A|. 
The number of spindown patches ~J\F S G S in the (s + 2)- 
dimcnsional search is 



NsG s 



s/2 



(a + 2)»/a 



(2.27) 



where Af s and G s are given in Eqs. (^15l) - (|2T9|) 1 and the 
prefactor on the right corrects for the s ky d imensions. 
The remaining product terms Y[ m Eq. ( 2.22| ) represent 
the increase in the size of the spindown space in order to 
include the orbital residuals. 



III. THRESHOLDS AND SENSITIVITIES 

The thresholds for a search are determined under the 
assumption that the detector noise is a stationary, Gaus- 
sian random process with zero mean and power spec- 
tral density S n (f). In the absence of a signal, the 
power P n (f) = 2|n(/)| 2 at each sampled frequency is 
exponentially distributed with probability density func- 
tion e~ Pn / Sn I ' S n . The statistic for stacked spectra is 
p = Pn(f)- The cumulative probability distribution 
function for p, in the absence of a signal, is 



rp/Sn 
CDF[p/S n ,N}= / 
Jo 



„N-1 



(N-iy. 



dr = 



l(N, P /S n ) 

(N-iy. 

(3.1) 



where "f(N, p/S n ) is an incomplete gamma function. 

A (candidate) detection occurs whenever p in some fre- 
quency bin exceeds a pre-specified threshold p c chosen so 
that the probability of a false trigger due to noise alone 
is small. There are / max AT Fourier bins in each spec- 
trum, and N p (AT, (j, max , N) spectra in the entire search. 
Therefore we assume that a search consists of N p f max AT 
independent trials of the statistic p, and compute the ex- 
pected number of false events F to be 
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F = / max AT N p (AT, /i max , N) (1 - CDF[p c / S n , N]) . 

(3.2) 

(In reality, there will be correlations between the statistic 
computed for different frequencies and different patches. 
Sinc e th is will reduce the number of independent trials, 
Eq. (^2) overestimates the number of false events. This 
is a small effect which should not change the overall sen- 
sitivity of a search by much. It is only in the case that 
the number of trials is initially small that one should be 
concerned with this effect; unfortunately, we operate in 
the other extreme.) If F <C 1, number of false events F 
is approximately equal to the probability that an event is 
caused by noise in the detector. Consequently, a = 1 — F 
can be thought of as our confidence of detection. In a 
non-hierarchical search, the threshold p c is set by speci- 
fying a and then inverting Eq. ([T^). 

Finally, how does the threshold p c affect the sensitivity 
of our search? We define a threshold amplitude /i t h to 
be the minimum dimensionless signal amplitude that we 
expect to register as a detection in the search 



^th = 



{ Pc /N - S n ) 



(F?(e,*,*))(l-</i>)AT 



(3.3) 




where (F^(0,$,^)) is the square of the detector re- 
sponse averaged over all possible source positions and 
orientations, and (/i) — /i m ax/3 is the expected mismatch 
of a signal which is randomly located within a patch. The 
sensitivity O of the search is then defined by 



(3.4) 



For any optimal search strategy, the goal of optimization 
will be to maximize the final sensitivity of the search, 
given limited computational power. 



IV. STACK-SLIDE SEARCH 

A stack-slide search is the simplest alternative to coher- 
ent searches we consider here. The main steps involved 
in the algorithm are shown in the flow chart of Fig. |l|. In 
this section we estimate the computational cost of each 
step, and determine the ultimate sensitivity of this tech- 
nique. 

The first step, before the search begins, is to specify 
the size of the parameter space to be searched (i.e. choose 
/max, t, and a region of the sky), the computational 
power P that will be available to do the data analysis, 
and an acceptable false alarm probability. From these, 
one can determine optimal values for the maximum mis- 
match /Lt max for a patch, the number of stacks N, and the 
length AT of each stack, using the optimization scheme 
discussed at the end of this section. For now, we treat 
these as free parameters. 



Coarse and fine grids are laid down on the param- 
eter space with N pc = N p (AT, /^ max , 1) and N p f = 
-/Vp(AT, /i max , N) points, respectively. The data-stream 
is low-pass filtered to the upper cutoff frequency f mgx , 
and broken into N stretches of length AT. Each of the 
steps above have negligible computational cost since they 
are done only once for the entire search. The subsequent 
steps, on the other hand, must be executed for each of 
the N pc correction points. 

Each stretch of data is re-sampled (at the Nyquist fre- 
quency 2/ max ) and simultaneously demodulated by stro- 
boscope sampling for a set of demodulation parameters 
selected from the coarse grid. The result is N demod- 
ulated time series, each one consisting of n = 2/ max AT 
samples. Since stroboscopic demodulation only shifts one 
in every few thousand data points (assuming a sampling 
rate at the detector of 16 384 Hz), the computational cost 
of the demodulation itself is negligible. 

Each stretch of data is then Fourier transformed using 
a fast Fourier transform (FFT) algorithm with a com- 
putational cost of 3nN log 2 (n) floating point operations. 
Power spectra are computed for each Fourier series, cost- 
ing 3 floating point operations per frequency bin, i.e., a 
total cost of 1.5nN floating point operations. 

For demodulation parameters in the coarse grid, the 
power of a matched signal will be confined to ~ 1 Fourier 
bin in each power spectrum, but not necessarily the same 
bin in different spectra. To insure that power from a 
signal is accumulated by summing the N spectra, we 
must apply N p f /N pc corrections within each coarse mesh 
patch. This can be achieved by the following steps. 

For the N spectra to be stacked, the frequency of a 
putative si gna l with initial frequency / max is computed 
using Eq. (1.1). Each spectrum is re- indexed so that the 
power from such a signal would be in the same frequency 
bin (we ignore the computational cost of this step), and 
the spectra are added [0.5n(iV— 1) floating point opera- 
tions] . We automatically account for corrections at other 
frequencies by applying the fine grid corrections in this 
way. It may be possible to reduce the computational cost 
of this portion of the search by noting, for example, that 
we over count the fine grid corrections for signals with 
frequency / max /2 by a factor of 2™ where n is the di- 
mension of the parameter space being explored. Since it 
is difficult to assess the feasibility of using this in a real 
search, we simply mention it so that it might be explored 
at the time of implementation. 

The resulting stacked spectrum is scanned for peaks 
which exceed the threshold p c . Since this has negligible 
computational cost, the number of floating point opera- 
tions required for the entire search is 

C = 3nNN pc [log 2 (n) + 0.5 + N pf (N - l)/(6NN pc )} . 

(4.1) 

If data analysis proceeds at the same rate as data ac- 
quisition, the computational power P required to com- 
plete a search is P = C/NAT floating-point operations 



11 




1 

0.8 
0.6 



w 
>> 
cd 
■a 



< 



10 11 10 12 10 13 10 14 10 15 10 16 
( a ) Computing Power (Flops) 

0.6 
0.5 

,0.4 
0.3 
0.2 
0.1 



III 1 1 III 1 III 


i mi inn i i 1 1 ii mi nil 

-'159 - 








/ 195 




/f92 




^192 


f X, 


6 


/1 92 








- /2A* 




mi i i 


i — 



10 



m 
cd 

[ 

< 



1 


1 1 II 1 1 1 Mil 


1 1 










^156 




- 


/236 


— 


: 




















/922 




- / 

m 


-^218 


1 ,1 



0.6 
0.5 
0.4 
0.3 
0.2 
0.1 

10 11 10 12 10 13 10 14 10 15 10 16 
( a ) Computing Power (Flops) 

0.6 
0.5 
,0.4 
0.3 
0.2 
0.1 



0.4 



0.2 ^ 

cd 

0.1 iH 

0.08^ 

0.06 

0.04 



II 1 1 III 1 1 III! 1 1 


11 1 1 1 1 11 1 1 1 1 1111 




/ 29 


f r 




/ 34 




: / 








- 

nil 1 1 





6 



in 
>> 
cd 

[ 

< 



10 7 10 8 10 9 10 10 10 11 10 12 
v°) Computing Power (Flops) 



10 11 1Q 12 1013 10 14 1Q 15 1Q 16 

(b) Computing Power (Flops) 



FIG. 7. The optimum stack length AT (thick solid line), 
number of stacks N (numbers along the solid line), and 
maximal projected mismatch /i max (dotted line) as functions 
of available computational power for directed, stack-slide 
searches. The two panels correspond to the following cases: 
(a) The situation encountered when searching for periodic 
sources having gravitational wave frequencies up to 1000Hz, 
with minimum spindown ages r m in = 40yr. (b) The equiv- 
alent results for gravitational wave frequencies up to 200Hz, 
with minimum spindown ages r m j n = 10 3 yr. Both cases as- 
sume a threshold which gives an overall statistical significance 
of 99% to a detection (although the results are insensitive to 
the precise value). The optimization was performed numer- 
ically using simulated annealing which accounts for some of 
the fluctuations in the observation times. 



per second (Flops). Equation (4.1) and the definition of 
n = 2/ max AT imply that the computational power is 

P = 6/ raax iV pc pog 2 (n) + 0.5 + N p f (JV - l)/(6JVJV pc )] . 

(4.2) 



The final sensitivity 0, defined in Eq. (3.4), of the 
search is determined once we know the function N p , the 



FIG. 8. Same as in Fig. (?], but for an all-sky search. 

frequency / max , the maximal mismatch /^ max , and the 
confidence level a = 1 — F. An optimized algorithm will 
maximize as a function of /x max , JV, and AT, subject 
to the constraints imposed by fixing the false alarm prob- 
ability F, and the computational power P. 

The results of the optimization procedure are given 
in Figures [7] and || for the fiducial classes of pulsar de- 
fined in Sec. |lB|. In each case we have set the proba- 
bility of a false alarm threshold at F = 0.01 (indicat- 
ing a 99% confidence that detector noise will not pro- 
duce an event above threshold), and have determined 
optimal values of /i max , JV, and AT for a range of val- 
ues of the available computational power. Figures 0(a) 
and (b) show the results for a directed search for old, 
slow (r m i n =l 000 yr, / max =200 Hz) and young, fast 
(T m i n =40 yr, /max=l 000 Hz) pulsars, respectively. Fig- 
ures ||(a) and (b) show the results for an all-sky search for 
the same two classes of source. The optimal sensitivities 
achieved by these searches are summarized in Fig. |^(a) 
of the Introduction. 
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V. HIERARCHICAL SEARCH: GENERAL 
REMARKS 

The basic hierarchical strategy involving a two pass 
search is represented schematically in Fig. 0. In the first 
pass, TVW stacks of data of length AT' 1 ) are demodu- 
lated on a coarse and fine mesh of correction points com- 
puted for some mismatch level /xP) , and then searched 
by stacked Fourier transforms. A threshold signal-to- 
noise level is chosen which will, in general, admit many 
false alarms. In the second stage, stacks of length 
are searched on a finer mesh of points computed 
at a mismatch level p^ , but only in the vicinity of those 
events which passed the first-stage threshold. The second 
stage will involve fewer correction points than the first, 
so the second-stage transforms can be made longer and 
more sensitive. The goal of optimization is to find some 



combination of AT«, ATP>, p^\ p (2 \ 7VP , 
which maximizes the final sensitivity for fixed computa- 
tional power P, and second pass false alarm probability 
F i2 \ 



A. Thresholds 



2F (D 



N$ AT™ 



T r(l-CDF\pW/Sg>\NW]) 



N W ATP 

= 2f max AT^N^(l - CDF[pP)/SP\ TVP)]) 



x(l-CDF[pW/sW,NW]) 



(5.2) 



where a is our desired confidence level for the overall 
search. 

The thresholds pP) and pP) cannot be assigned inde- 
pendently; rather, they should be chosen so that any true 
signal buried in the noise that would exceed (in expec- 
tation value) the second-stage threshold will have passed 
the first-stage threshold. In other words, it serves no 
purpose to set pP) any lower than the weakest signal 
which would have passed pP). A signal which is ex- 
pected to pass the second-stage threshold exactly has an 
and™ amplitude \U^\ 2 = p« - N^S^. We define the false 



dismissal probability D to be the probability that such 
a signal will be falsely rejected in the first pass. Since 
the spectral power of a true signal increases with NAT, 
the signal seen in the first pass has amplitude |/z,P^| 2 = 
\h( 2 )\ 2 (N^AT^)/(N^AT^), and the thresholds sat- 
isfy the relation 



In the first pass of a hierarchical search, each 
of = / max ATP) frequency bins in N$ = 

iVp(ATP),^P),JVP)) stacked power spectra will be 
scanned for threshold crossing events. If (as we assume) 
all of these trials are statistically independent, the num- 
ber of false events above the threshold pP' will be 

F 0.) = Jvg>jvW(l - CDF[pP)/5W,ArP)]) . (5.1) 

We assume that the number of false events will signif- 
icantly exceed the number of true signals in this pass, 
consequently the number of events to be analyzed in the 
second pass will be 

The second stage uses a coarse grid with N p t 
N p (AT (2 \p (2 \ 1) points, and a fine grid with N f 
N p (AT (2 \p (2 \nW) points. On average each false 
alarm will require Npc /N^ 1 ) coarse grid points, and 

N^j /N^j fine grid points in the second stage. (When 
a second-pass mesh is coarser than the first pass's pa- 
rameter determination, the corresponding ratio should 
be taken as unity.) Furthermore, since the first stage will 
identify the candidate signal's frequency to within ~ 2 
frequency bins, the second-stage search should be over 
the 2AT( 2 )/ATP) second-stage frequency bins which lie 
in this frequency range. Once again, we assume the noise 
in all frequency bins (and over all grid points) is inde- 
pendent, so the number of false events which exceed the 
threshold pP) in the second stage is 



(2) 
>c 
(2) 



1 



D = CDF 



= CDF 



c-(l) 
On 



-,atP) 



c(l) 
On 



pP) \ ^iV«ATP) , 

0(2) I cP) ATP)ATP)' 

o n / On 



Now, for any choice of AT^\ AT {2 \ etc., the thresholds 
pP) and pP) are completely constrained by our choices 
of final confidence level a and false dismissal probability 
D. The false dismissal probability is fixed at D — 0.01 in 
our optimization; this is an acceptably low level, meaning 
that only one signal in a hundred is expected to be lost 
in this type of search. 



B. Computational costs 

The computational cost CP-* of the first stage of the 
search follows the same formula as for a simple non- 
hierarchical search, that is 

C (1) = 6/ max iVP)ATP)iVP) [log 2 (2/ max ATP)) 

+0.5 + ATP)(jVP) - 1)/(6JVP)JVP))] . (5.4) 

For each of the F^ first-stage triggers, the second 
stage requires Npc /N^j) (minimum 1) coarse grid cor- 
rections (each involving N^ 2 ' FFT's of length ATP)), 

along with N^y /N^ 1 ) (minimum 1) frequency shifts and 
spectrum additions. Each of the coarse grid corrections 
requires the usual 2/ max iVP) ATP) [3 log 2 (2/ max ATP)) + 
0.5] floating-point operations. The incoherent frequency 
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shifts and spectrum additions require only 2(iVd) — 
V) AT^ 2 ' / AT^ 1 ' floating point operations since the fre- 
quency correction and power summation need only be 
applied over a bandwidth of ~ 2 first-pass frequency bins. 
The total cost of the second pass is therefore: 
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FIG. 9. The optimum stack length AT (thick solid line) 
and number of stacks TV (numbers along the solid line) for 
the first and second stages of directed, hierarchical searches. 
For numerical convenience the maximal projected mismatch 
was chosen in advance to be /i max = 0.3, and the number 
of stacks in the first and second stage are the same. The 
two panels correspond to the following cases: (a) The situa- 
tion encountered when searching for periodic sources having 
gravitational wave frequencies up to 1000Hz, with minimum 
spindown ages r m i n = 40yr. (b) The equivalent results for 
gravitational wave frequencies up to 200Hz, with minimum 
spindown ages r m j n = 10 3 yr. Both cases assume a thresh- 
old which gives an overall statistical significance of 99% to a 
detection (although the results are insensitive to the precise 
value) . 



C (2 



N 



(i) 
pf 



3/max[l0g 2 (2/ max AT( 2 )) 



-0.51 



nwn$atw 



(5.5) 



We require that data analysis proceed at the rate of 
data acquisition. Since the amount of data used in the 
second-stage of the search will generally be greater than 
that used in the first, we require that the analysis be com- 
pleted in N^AT^ seconds. Thus the computational 
power is given by 



P= (C^ +C^)/N^AT^ 



(5.6) 



Our final sensitivity is given by Eq. (3.4), using 
the observation time, mismatch level, and threshold of 
the second stage of the search. Optimization then con- 
sists of maximizing this function over the six parameters 
AT«, AT( 2 ), ^(D, /id), /V«, and given constraint 



Eq. (5.6) for specified a, D, and P. 



VI. HIERARCHICAL SEARCH WITH STACKING 

It turns out that the optimization described in the pre- 
vious section is only weakly sensitive to the parameters 
/«d) and /td) j that is, even if we choose values for /«d) and 
/id) quite different from the optimal ones, we can recover 
nearly all of the sensitivity by adjusting the other param- 
eters for the same computational power P. In particular, 



if we arbitrarily fix \x 



(i) 



(2) 



0.3 and re-optimize, we 



obtain sensitivities within 20% of the optimal. 

This becomes very useful when we consider the gener- 
alized two-stage hierarchical search with stacking. Nor- 
mally this would involve optimizing over six variables 
(/id)<d), /vO)<( 2 ), and ATd)>d)) with one constraint on 
P. However, by assuming that we can continue to set 
/id) = /id) = 0.3 with minimal loss of sensitivity, we can 
reduce our degrees of freedom back down to four minus 
one constraint. 

The results of this optimization for our four canoni- 
cal example searches are given in Figs. ^ and [l(]. We 
have again chosen a final confidence level a = 0.99 
and a false dismissal probability of D = 0.01. Fig- 
ures ||(a) and (b) show the results for a directed search for 
old, slow (r m in=l 000 yr, / max =200 Hz) and young, fast 
(T m in=40 yr, / max =l 000 Hz) pulsars, respectively. Fig- 
ures [n](a) and (b) show the results for an all-sky search 
for the same two classes of source. The optimal sen- 
sitivities achieved by these searches are summarized in 
Fig. ||(b) in the Introduction. 
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FIG. 10. Same as in Fig. Bl but for an all-sky search. 



VII. SPECIALIZED SEARCHES 

The strongest sources of continuous gravitational 
waves are likely to be the most difficult to detect since 
the frequency of the waves will be changing significantly 
as the source radiates angular momentum. As we have 
seen in the previous sections, an all sky search for these 
sources is unlikely to achieve the desired sensitivity with 
available computational resources. To reach better sensi- 
tivity levels, it will be useful to consider targeted searches 
for specific types of source. In this section, we consider 
three such searches: (i) neutron stars in the galactic core 
as an example of a limited area sky survey, (ii) new- 
born neutron stars triggered on optically observed extra- 
galactic supernovae, and (iii) low mass X-ray binary sys- 
tems such as Sco X-l. 



A. Galactic core pulsars 

Area surveys of the sky will certainly begin with the 
region most likely to hold a large number of nearby 
sources. Based on population models of radio pulsars 
in our Galaxy |30|, there should be many rapidly rotat- 
ing neutron stars in the galactic bulge. As an exam- 
ple of a limited area search, we therefore consider the 
optimal strategy for searching an area of 0.004 steradi- 
ans about the galactic core, for sources with frequencies 
/ < 500 Hz and spindown ages r > 100 yr. The choice 
of a 0.004 steradian search is arbitrary; it includes the 
entire molecular cloud complex at the core of the galaxy 
(~ 300 pc radius at a distance of ~ 8.5 kpc). 

It is easy to include a correction factor, to allow for 
this limited area, in our calculation of the number of 
patches by reducing the ranges of the integral over V in 
Eq. (2.14). Given the approxi mation s in Sec. [IB, this 
amounts to reducing N p in Eq. ( |2.22j ) by 

where the multiplicative factor 0.97 is the correction for 
the difference in functional form between the mismatch 
metric and the angular area metric dSl 2 — ; " 
the direction of the galactic center (i.e. 
tion). 

The optimal choices of 2VW>( 2 ) and AT^'^ for a hier- 
archical stacked search are shown in Fig. [ll] as a function 
of available computing power; the relative sensitivity of 
this search is shown in Fig. ^| of the Introduction. 
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FIG. 11. The optimum stack length AT (thick solid line) 
and number of stacks N (numbers along the solid line) for the 
first and second stages of an hierarchical search for pulsars 
located in a sky region of 0.004 steradians about the Galactic 
center, with r > 100 yr and / < 500 Hz. For numerical 
convenience the maximal projected mismatch was chosen in 
advance to be /x max = 0.3, and the number of stacks in the first 
and second stage are the same. A threshold is chosen which 
gives an overall statistical significance of 99% to a detection 
(although the results are insensitive to the precise value). 
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We note from Eq. (3.6) of Paper I that gravitational 
waves from rapidly rotating neutron stars might be ex- 
pected to have a characteristic amplitude of 



2.3 x 10 



-25 



I, 



8.5kpc / / 



10-5 10 45 gcm 2 



500Hz 



(7.2) 



where e = (J IX - I yy ) / 'I zz is the non-axisymmetric strain, 
ijj- is the moment of inertia tensor, r is the distance to 
the source, / is the gravitational wave frequency, and h c 
has been averaged over the detector responses to vari- 
ous source inclinations ]27j . Theoretical estimates of the 
strength of the crystalline neutron star crust suggest that 
it can support static deformations of up to £ ~ 10~ 5 , 
though most neutron stars probably support smaller de- 
formations. From Figs. || and [g, we see that 1 Tflops of 
computing power should allow us to detect pulsars with 
strains as small as e ~ 5 x 10 -6 at 8.5kpc using enhanced 
LIGO detectors. 



B. Newborn neutron stars 

Several recent papers ]2p|-|22]| have indicated that 
newly-formed fast-spinning neutron stars may be copious 
emitters of gravitational radiation. If the newborn neu- 
tron star is rotating sufficiently fast, its r- modes (axial- 
vector current oscillations whose restoring force is the 
Coriolis force) are unstable to gravitational radiation re- 
action. As the star cools, viscous interactions eventu- 
ally damp the modes in isolated neutron stars. Numer- 
ical studies |23[ indicate that neutron stars which are 
born with rotational frequencies above several hundred 
Hz will radiate away most of their angular momentum in 
the form of gravitational waves during their first year of 
life. Estimates of the viscous time scales, and the super- 
fluid transition temperature, suggest that the r-modcs 
are stabilized when the star cools below ~ 10 9 K and 
are rotating at ~ 100-200 Hz. During the evolutionary 
phase when most of the angular momentum is lost, the 
amplitude and spindown time scale are expected to be 



580s ZlkHzV 



6i 



(7.4) 



These estimates are based on Eqs. (4.9) and (5.13) in 
Rcf. p3| . (We note that the "characteristic amplitude" 
used in Ref. Q is appropriate to estimate the strength 
of burst sources, and is different from our h c .) Here n is 
a dimensionless constant of order unity; it parameterizes 
our ignorance of the non-linear evolution of the r-mode 
instability. The distance to the neutron star is r, and 
t is the actual age of the star. Figure || shows h c as a 



function of frequency with k = 1 at distances r = 2 Mpc 
and 20 Mpc. 

Sources outside our Galaxy are potentially detectable 
due to the high gravitational-luminosity of a newborn 
neutron star with an active r-mode instability. Never- 
theless, it is a significant challenge to develop a feasi- 
ble search strategy for these signals since the frequency 
evolves on such short time scales (compared to those 
considered above). One approach is to perform directed 
searches on optically observed supernova explosions. Al- 
though some supernovae may not be optically visible, and 
this instability may not operate in all newborn neutron 
stars, the computational benefits of targeting supernovae 
are substantial (if not essential). Based on the estimates 
in Ref. , most of the signal to noise is accumulated 
during the final stages of spindown. With limited com- 
putational resources, it seems best to limit the directed 
searches to frequencies < 200 Hz, when the spindown 
time scale is ~ 1 yr. Figure |l2| shows the optimal search 
criteria in a hierarchical stacked search for neutron stars 
aged two months or older; the upper frequency cutoff is 
/max = 200 Hz and the minimum spindown timescale is 
Tmin = 1 yr. The sensitivities achievable in a search are 
shown in Fig. |^ of the Introduction. 
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FIG. 12. The optimum stack length AT (thick solid line) 
and number of stacks N (numbers along the solid line) for the 
first and second stages of a hierarchical search for newborn 
neutron stars spinning down due to an active r-mode insta- 
bility. We assume that a supernova has been identified and 
accurately located on the sky, so this is a directed search for 
an object with r m i n = 1 yr, and / < 200 Hz. For numeri- 
cal convenience the maximal projected mismatch was chosen 
in advance to be fi m&K = 0.3, and the number of stacks in 
the first and second stage are the same. A threshold is cho- 
sen which gives an overall statistical significance of 99% to a 
detection (although the results are insensitive to the precise 
value) . 

Figure [| shows that 1 Tflops of computing power will 
not suffice to detect newborn neutron stars as far away as 
the Virgo cluster (« 20 Mpc); however, such sources will 
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be marginally detectable within ~ 8 Mpc by enhanced 
LIGO detectors. The NBG catalog §lj lists 165 galax- 
ies within this distance (assuming a Hubble expansion 
of 75 km/s/Mpc, retarded by the Virgo cluster). From 
the Hubble types and luminosities of these galaxies, and 
the supernova event rates in |^| , we estimate a total su- 
pernova rate of ~ 0.6 per year in this volume, of which 
~ 10% would be of type la, ~ 20% of type lb or Ic, and 
~ 70% of type II. (We note that the total rate is consis- 
tent with values given in Ref. |p3fl .) At present, it is not 
known what fraction of these will produce neutron stars 
with unstable r-modes. 



C. X-ray binaries 

A low-mass x-ray binary (LMXB) is a neutron star 
orbiting around a stellar companion from which it ac- 
cretes matter. The accretion process deposits both en- 
ergy and angular momentum onto the neutron star. The 
energy is radiated away as x-rays, while the angular mo- 
mentum spins the star up. Bildsten pij has suggested 
that the accretion could create non-axisymmetric tem- 
perature gradients in the star, resulting in a substantial 
mass quadrupole and gravitational wave emission. The 
star spins up until the gravitational waves are strong 
enough to radiate away the angular momentum at the 
same rate as it is accreting; according to Bildsten's esti- 
mates the equilibrium occurs at a gravitational- wave fre- 
quency ~ 500 Hz. The characteristic gravitational-wave 
amplitudes from these sources would be 



h c > 4 x 10" 



R 



10km 
F 



3/4 



M 



10 _8 erg cm _2 s" 



1AM® 

1/2 



-1/4 



( f 



V 600 Hz 



-1/2 



, (7.5) 



where R and M are the radius and mass of the neutron 
star, and F is the observed x-ray flux at the Earth. 

The amplitude of the gravitational waves from these 
sources make them excellent candidates for targeted 
searches. If the source is an x-ray binary pulsar — an ac- 
creting neutron star whose rotation is observable in radio 
waves — then one can apply the exact phase correction 
deduced from the radio timing data to optimally detect 
the gravitational waves. (In this process, one must as- 
sume a relationship between the gravitational-wave and 
radio pulsation frequencies.) Unfortunately, radio pulsa- 
tions have not been detected from the rapidly rotating 
neutron stars in all LMXB's (i.e. neutron stars which 
rotate hundreds of times a second) . In the absence of di- 
rect radio observations, estimates of the neutron-star ro- 
tation rates are obtained from high-frequency periodic, 
or quasi-periodic, oscillations in the x-ray output dur- 
ing Type I x-ray bursts. (See Ref. j34|] for a summary.) 
But this does not provide precise timing data for a co- 
herent phase correction. To detect gravitational waves 



from these sources, one must search over the parame- 
ter space of Doppler modulations due to the neutron- 
star orbit around its companion, and fluctuations in the 
gravitational-wave frequency due to variable accretion 
rates. The Doppler effects of the gravitational-wave de- 
tector's motion can be computed exactly, because the sky 
position of the source is known. 

In most cases, the orbital period of an x-ray binary can 
be deduced from periodicity in its x-ray or optical light 
curve. In some cases, the radial component of the orbital 
velocity can be computed by observing an optical emis- 
sion line from the accretion disk, as was done with the 
bright x-ray binary Sco X-l Such observations do 

not determine the phase modulation of the gravitational- 
wave signal with sufficient precision to make the search 
trivial; however, they do substantially constrain the pa- 
rameter space of modulations. 

In this subsection, we consider a directed search for an 
x-ray binary in which the orbital parameters are known 
up to an uncertainty Sv in the radial velocity v r of the 
neutron star, and an uncertainty S(f> in the orbital phase. 
It is assumed that long-term photometric observations 
of the source can give the orbital period P to sufficient 
precision that we need not search over it explicitly. We 
therefore parameterize the phase modulations as follows: 

M; A) = 2tt/o (t + — cos 2nt/P + — sin 2nt/P ) , 

\ 2lTC 2lTC ) 

(7.6) 

where A = (/o,Ui,i>2) are our search parameters, the 
gravitational-wave frequency fo is constrained to be < 
/max, and the pair (v\,V2) is constrained to lie within an 
annular arc of radius v r , width Sv, and arc angle 5<j>. 

Applying the formalism developed in section [H] to this 
problem gives essentially the same result as for a sky 
search over Earth-rotation-induced Doppler modulations 
if one converts time units by the ratio P/day. In the 
case of the Earth's rotation, a search over sky positions 
h corresponds to a search over an area Tiv^ ot cos 2 (A) in 
the equatorial components of the source's velocity rela- 
tive to the detector, whereas in the case of a binary orbit, 
the search is over a coordinate area v r SvS(f>. So we can 



simply multiply equation (2.22) by the ratio of these co- 
ordinate areas to obtain the number of grid points N p in 
the parameter space: 



( /max P) 2 ^ ( ^ 2+i3 _ 2+c _ 2) -l/2 ) (7?) 



A = 0.5 , 

B = 1.6(AT/P) 2 , 
C = 6A{AT/P) 5 N 3 



(7.8) 
(7.9) 
(7.10) 



Accounting for the intrinsic phase variations of the 
spinning neutron star itself is problematic. Baykal and 
Ogelman p6| showed empirically that x-ray pulsar fre- 
quencies could be well-modeled as a random walk, plus 
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possibly a secular spin up for rapidly-accreting sys- 
tems. Over a time t, the angular rotation rate would 
undergo excursions up to Afl = \/Si, where S ~ 
10 _17 rad 2 s _3 (L x /10 37 ergs _1 ) and L x is the x-ray lu- 
minosity of the source. For the sources of interest to 
us, accretion proceeds at or near the Eddington rate 
(L x ~ 3 x 10 38 erg/s), and the gravitational- wave fre- 
quency is / = f2/7T, so we expect frequency drifts ~ 
2 x 10~ 6 Hz^/t/days. If we require that the frequency 
drifts by less than one Fourier bin A/ < 1/t during a 
coherent observation, the observing time t must satisfy 

t< 3.2 days. (7.11) 
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FIG. 13. The optimum stack length AT (thick solid line) 
and number of stacks N (numbers along the solid line) for the 
first and second stages of a hierarchical search for sources of 
continuous gravitational waves which are in binary systems. 
We assume the orbit is characterized by two orthogonal ve- 
locity parameters which are known to within a total error of 
17(km/s) 2 , and that the frequency / < 500 Hz is experienc- 
ing a random walk typical of Eddington-rate accretion. For 
numerical convenience the maximal projected mismatch was 
chosen in advance to be /i max = 0.3. A threshold is cho- 
sen which gives an overall statistical significance of 99% to a 
detection (although the results are insensitive to the precise 
value) . 

This type of random walk cannot be modeled as a low- 
order polynomial in time. Nevertheless, the stack-slide 
technique is well suited to search for these sources since 
corrections for the stochastic changes in frequency can be 
applied by shifting the stacks by +1, 0, or —1 frequency 
bins when required. In a search using N stacks, each of 
length AT, this kind of correction would be applied after 
a time t such that 

2 x HT 6 HzVt/days = 1/AT . (7.12) 

The number of times that these corrections must be ap- 
plied is then NAT/t, and the number of distinct fre- 
quency evolutions traced out in this procedure is 3 JVAT /*. 



(Monte-Carlo simulations of stacked FFTs of signals un- 
dergoing random walks in frequency have shown that one 
can increase t by up to a factor of 4, i.e. allowing drifts 
of up to ±2 frequency bins, while incurring only ~ 20% 
losses in the final summed power.) We have not yet stud- 
ied in detail how this combines with the mismatches gen- 
erated from other demodulations, or how to search over 
all demodulations together in an optimal way. For now 
we assume a factor of 3° o lA, ( AT / da y) extra points in our 
search mesh for mismatches of ^ max — 0.3. 

As an example, we consider a search for gravitational 
waves from the neutron star in Sco X-l. This system has 
an orbital period P = 0.787313 ± 0.000001 days |37), a 
radial orbital velocity amplitude of v r — 58.2 ±3.0 km/s, 
an orbital phase known to ±0.10 radians |3f|, and an 
inferred gravitational- wave frequency / ma x ~ 500 Hz [ p4[ . 
We note that the uncertainty in P is basically negligible 
over the < 3.2 day coherent integrations expected. The 
remaining uncertainties give: 

N p = —3 omN ( AT / d ^ 3 (A- 2 + B- 2 +C- 2 y 1/2 , (7 

A = 0.5 , (7 
B = 2.6(AT/day) 2 , (7 
C = 21Ar 3 (AT/day) 5 . (7 

where it is understood that AT < 3.2 days, in order for 
the random-walk stack-slide corrections to achieve max- 
imum sensitivity. 

Figure [ll| shows the optimal search criteria for a hi- 
erarchical stacked search for the Sco X-l pulsar under 
these assumptions. The sensitivities achievable in such 
a search are shown in Figure ^| of the Introduction. We 
see that 1 Tflops of computing power may be sufficient 
to detect this source using enhanced LIGO detectors if 
it is radiating most of the accreting angular momentum 
as gravitational waves. The sensitivity to these sources 
might be enhanced by a factor ~ 5-10 if the interferom- 
eter is operated in a signal-recycled, narrow-band config- 
uration during the search. 

VIII. FUTURE DIRECTIONS 

We have presented in this manuscript the rudiments of 
a search algorithm for sources of continuous gravitational 
waves. The next step is to implement some variant of 
these schemes for a simple test search; a good starting 
point would be a directed, acceleration (a single spindown 
parameter) search of a nearby supernova remnant. This 
will be sufficient to test every stage of the technique, and 
to assess the accuracy of the computational estimates 
presented here. It will also allow direct comparison with 
other approaches such as the line-tracking method that 
is being explored by Papa Jl6| . 

Further theoretical work is required to determine the 
parameter space which should be searched, especially in 
the case of active r-mode instabilities and radiating neu- 
tron stars in LMXB's. Finally, it would be worthwhile to 



18 



consider in detail what advantage, if any, can be gained 
by using data from multiple interferometers at the initial 
detection stage of a search for continuous gravitational 
waves. 
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APPENDIX A: PATCH NUMBER FORMULAE 



The approximate formulae given in Eqs. (pl7j)- (pls|) are valid when N 4. General expressions for the G s can 
be derived by setting x = in Eqs. (Q and (|l 6|), and using Eqs. @, (fEj), ( pU0| ), ( pH] ), ( ]2~13| ), and ( pTT^ ) : 



G 2 (iV) 



V1757V 6 - 840JV 4 + HOOiV 2 - 432 



(Al) 
(A2) 



180V105 

G 3 (N) = = V36757V 12 - 58800JV 10 + 3631607V 8 - 10533607V 6 + 14843367V 4 - 9878407V 2 + 248832 . (A3) 



75600V105 

In Eq. ( |2. 14 ) we approximate the metric as having constant determinant and evaluate it at the point of zero spindown; 
this introduce s sm all erro rs of order ([NAT]/t) 2 . 

Equations (2.22 )-( 2.26| ) for the number of sky patches ignoring spindown and orbital motions provide an empirical 
fit to a numerical evaluation of the metric determinant. The determinant was found to have an approximate func tiona l 
depe nden ce ^/pj^j | ~ | sin2#| with corrections of order v/c ~ 10~ 4 . Assuming this dependence to be exact, Eqs. (2.12) 
and ( 2.14 ) give 



Ms = 



2/tt 



4aw x /(s + 2) 



\ 



det 



.90' 



9oi9oj 



9oo 



« = ir/2 



(A4) 



Here g a p is the mismatch metric, defined in Eqs. (|2.9| )-( ^1~1 ), computed using only the Earth-rotation- induced Doppler 
modulation. Since the Earth's rotation is a simple circular moti on, an d since we are evaluating the metric at a single 
point in parameter space, we can carry out the integrals in Eq. ( 2.11 ) analytically, to obtain 



1 N 

Ja(3 — b afi - — a ka 



a-k/3 



k=l 



where 



am 

Ofc2 

boo 
hi 
bi2 

&01 
&02 



1 

~ 2 

1 

12$ 



v/c 



{sin(fcA$) - sin([fc - 1]A$)} 



{sin(fcA$) - sinflfc - 1]A$)} 
{cos([fc - 1]A$) - cos(/cA$)} 



24 + - + 6 - $ + 4$ 



24- cos $ + 24-$ sin $ + 3 

c c 



(7 



sin 2$ 



/L("/c) 



4$ 

/L(«/c) 2 

4$ 

7 /max U 

010 

&20 



(2$ + sin 2$) 

(2$ - sin 2$) 

<^ -4 + 2-$ + 4 cos $ + 4$ sin $ + - sin 2$ } 
4$ I c c J 

•^ m ^ C | -2$ cos $ + sin $ (2 + - sin j 



6 12 = 6 21 = fh»M^? sin 2 $ . 



(A5) 

(A6) 

(A7) 
(A8) 
(A9) 

(A10) 
(All) 
(A12) 
(A13) 
(A14) 



Here $ = 27rTf,/(l day) is the total angle over which the Earth rotates during the observation, A$ = Q/N is the 
angle rotated during each stretch of the data, and v is the maximum radial velocity relative to the detector at latitude 
A = 45° of a point at a polar angle 9 — 7r/4 on the sky, that is 



27T_R eart h cos A sin 9 
lday 



(A15) 
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APPENDIX B: RESAMPLING ERROR 

In this paper, we have assumed that coherent phase 
corrections are achieved through stroboscopic resam- 
pling: a demodulated time coordinate tf, [t] is constructed, 
and the data stream h(t) is sampled at equal intervals 
in tb at the Nyquist rate for the highest frequency sig- 
nal present, /Nyquist = 2/ max . However, since the data 
stream is initially sampled at some finite rate f s — 
R /Nyquist (where R is the over sampling factor), this can 
introduce errors: in general, there will not be a data point 
exactly at a given value of tb, so the nearest (in time) da- 
tum must be substituted. Consequently, there will be 
residual phase errors A<i>(i) g [—ir/2R,ir/2R) caused by 
rounding to the nearest datum even if one chooses a phase 
model whose frequency and modulation parameters ex- 
actly match the signal. The phase of the resampled signal 
drifts until the timing error is > l/2/ s , at which point 
one corrects the phase by sampling an adjacent datum, 
which shifts in time by l// s . These residual phase errors 
reduce the Fourier amplitude of the signal by a fraction 



F 



1 N 

fe=i 



(Bl) 



where N — /Nyquist AT is the total number of points in 
the resampled data stream. 

The uncorrected signal will in general drift by many 
radians in phase, which is the reason why we must ap- 
ply phase corrections in the first place. This means that 
A$(t) will sweep through the range [— ir/2R, ir/2R) many 
times over the course of the observation. So, regardless 
of the precise form of the phase evolution, we expect 
A $(fc/ /Nyquist) to be essentially evenly distributed over 
this interval. Thus, replacing the sum with an expecta- 
tion integral, we have 



F 



TT/2R 



n/2R 



sin(7r/2i?) 
n/2R 



(B2) 



The fractional losses in amplitude 1 — F for a few values 
of the oversampling factor R are given in Table |. 

TABLE I. The percentage reduction (1 — F) in amplitude 
of a signal as a function of the oversampling factor R. The 
LIGO interferometers will collect data at f s = 16 384Hz, so 
that the data will be oversampled by R > 4 compared to 
the maximum gravitational wave frequency that we expect 
on physical grounds. In fact, it seems more likely that R ~ 8 
for real signals. 



R= 


2 


3 


4 


5 


6 


7 


8 


1-F= 


10.0% 


4.5% 


2.6% 


1.6% 


1.1% 


0.8% 


0.6% 



The highest gravitational- wave frequencies we consider 
are 1000 Hz, requiring a resampling rate of /Nyquist = 
2000 Hz. Since LIGO will acquire data at a rate of 



16 384Hz, corresponding to an oversampling factor of 
R > 8, we have a maximum signal loss due to resampling 
of 1 — F = 0.6%. Resampling errors will increase if the 
number of data samples is reduced by some factor before 
phase-correcting. 
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